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ABSTRACT 


A numerical solution for the KdV solitons has been 
obtained using Galerkin scheme. The basis functions used 
are cubic splines. The scheme gives a system of ordinary 
differential equations in the time variable. The solution 
in time is evolved using IMSL library subroutines. 

The solution obtained is compared with the various 
finite difference and finite element schemes used previously. 
The errors are found to be significantly less as compared 
with the finite difference schemes reported in the literature. 
The Galerkin scheme used here is also found to give better 
results . than the Petrov-Galerkin and the Modified Petrov- 
Galerkin schemes reported in recent literature. 

The results obtained here also appear to be generally 
more favorable than the dissipative Galerkin schemes based on 
cubic splines and quintic boundary basis functions used 
previously. 
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CHAPTER 1 
INTRODUCTION 


1 .1 Forewerd 

It is well known that the initial and boundary value 
problems associated with nonlinear partial differential 
equations are very difficult to handle in a general way. 

Some speeific problems have been tackled from time to time 
by methods specifically suited to the individual problems. 

Over the past two decades, a remarkable development 
in our understanding of certain classes of nonlinear partial 
differential equations known as (nonlinear) evolution 
equations has taken place. A fundamental relationship 
between some of these nonlinear equations and ordinary 
linear differential equations of Sturm-Liouville type has 
been demonstrated through a variety of Spectral Transform, 
Even at our present level of understanding, the basis for 
such a fundamental relationship remains somevtiat obscure [l3. 

The main reasons of interest in these nonlinear 
evolution equations are two fold: 

1 . These arise in a natural way in a large number 


of physical problems. 
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2, They possess special type, of solutions \Atiich 
can be of great practical use. 

These special solutions take the form of localized 
travelling disturbances, or pulse, with peculiar properties, 

1, Their speed, (effective) width and amplitude 
are inter-related unlike the solutions of linear partial 
differential equations which have travelling solutions, 

2, Their shape does not get distorted as they travel 
through the medium, 

3, They retain their shape, amplitude and speed (i.e, 
complete identity) even after interactions among themselves, 
and suffer only a phase shift upon such interactions. 

Due to above mentioned particle-like properties, 
such localized travelling disturbances have come to be known 
as Solitons. 

Our interest in the present work is in the numerical 
simulation studies of single-sol it on solutions of one such 
evolution equation known as Korteweg-deVries (KdV) equation. 
The KdV equation arises in a number of hydrodynamic problems 
including problems relevant to plasma physics. 

In the remainder of this chapter, we introduce the KdV 
equation and briefly discuss some physical situations leading 
to the KdV equation. This is followed by a brief review of 
the literature and an outline of the present work. 
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1 .2 The Korteweq- deVries (KdV) Equation and Its Transformations 
The KdV equation may be written as; 


A 


iu 

3t 


Bu 


o u 

3' X 


+ C 




= 0 , 


(1.1) 


where u = u(x,t), x and t are independent variables, and 
A, B and C are all real, nonzero constants. Upon division 
by A, one can wTite Eq.(l,l) as; 



where a = B/A and b = C/A. It should be pointed out that 
suitable transformations of u, x and t allow the introduction 
of arbitrary constants multiplying the three terms in Eq,(l.l) 
or ( 1 ,2) , For example, a rescaling of Eq,(1,2) with 


will yield: 


u -♦ u/ a 


Ut + uuj^ + e u^^^ = 0 , 


( 1 .3) 


where e= b. In writing Eq,( 1 ,3) we have used a standard 
notation for the partial derivatives, viz,; 


Uj 


8u 


“x 2 ^ and 


» u 

r? 

0 A 


A rescaling with 


u 


1 /3 

^ X b , and 

„ o-l hl/S 
-* u a b 


will give coefficients of unity in front of each terms, i,e, , 



4 


Uj 


+ uu 


+ u 


XXX 


(1.4) 


1/3 1 *“1 1/3 

Instead, a transformation x xb ' and u ~ u a b 

in Eq,( 1 ,2) or a rescaling with u ~ u in Eg, (1.4) will 

yield another useful form of the KdV equations 


“t - + “xitx = ^’-5) 

As will be seen in Section 2.2, the form (1.5) above is most 
convenient for a discussion of the Spectral Transform Method. 

It may be mentioned that some of these transformations 
may not be necessarily unique. For example, in Appendix A, 
we give two different transforms which change Eq,(l.l) to 
Eq,( 1 ,14) . 

It is perhaps worthwhile to point out that if s < 0 in 
Eq.(1.3), a transformation u -u and x -* -x will lead to 
a positive coefficient of the term u^^^, 

X X. X 

In the present work we would be interested in the 
solution of the KdV equation for localized initial data 
u(x,0), where u(x,0) and its derivatives with respect to x 
vanish sufficiently rapidly as j x 1 -» 00 . It is also of 
interest to note that without the nonlinear term, uu , Eq.(1.3) 

A 

describes linear dispersive waves, whereas, without the 
dispersive term, u , the same equation describes the well 
known shock phenomena in fluid dynamics [2,3)1. 
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1 ,3 Physical Systems Leading to the KdV Equation 

As has been remarked earlier a number of physical 
systems, mostly related to hydrodynamics, lead to the KdV 
equation. Unfortunately, somewhat tedious perturbation 
calculations are needed before the equation of interest, 
such as the KdV equation, is finally obtained. We briefly 
mention here two examples, 

1 .3,1 Shallow Liquid Waves 

We consider the motion of plane disturbance on the 
surface of a liquid in a channel of constant depth. 

Let, 

h = the constant depth of the liquid in the channel 
a = the amplitude of the disturbance, and 
1 = the effective width of the disturbance along 

the channel. 

If one assumes that 

■| < < 1 , and 

< < 1 , 

a perturbation development of the resulting fluid motion can 
be obtained. If the fluid motion is assumed to be nonviscous, 
incompressible, irrotational, and bounded below by a hard 
horizontal bed, then the following equation for the resulting 
motions is obtained [1»4]: 



6 


2 K 


° TTH + i c crri 

X 2 0 XXX 


( 1 . 6 ) 


where 

T) ^ h (5<»t), the boundary of the free surface. 



g = acceleration due to gravity, 
2 

^ ^ T 

o - “3" » 

T = surface tension, and 

p = density of the fluid. 


2 

If the surface tension is negligible, i.e, T << h pg, the 

1 1 2. 

quantity a in Eq,(1 .6) simplifies to ^ • [1] « 

It is now straightforward to see that the Galilean 
tr ansf ormat ion 


X -* X + c^t 

transforms Eq.(1,6) to the KdV Eq.(1,2), 

1 .3.2 Ion Plasma VVaves and the KdV Equation 

In this section we will show that nonlinear ion- 
acoustic waves travelling with a speed close to the linear 
acoustic speed in a cold plasma (T^ = 0) obey the KdV 
equation and thus possess the soliton characteristics 
consider the simple case of one-dimensional disturbance. 
The fluid equations of motion and continuity are [ 5,6,7 3 
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3 V. 3 V. 

^ + Vi ^ 


e fjP 
m 3x 


(1.7) 


3ni 3 

cnr” **’ "Sjc = 0 (1 ,8) 

where 

Vi = speed of the ions, 
e = magnitude of electron charge, 
m = electron mass, 

9 *= electric potential, 

ni = ion density, and, 

X and t denote the space and time variables respecti- 
vely, If we neglect the electron inertia, the Poisson's 
equation takes the form [5] 


,2 e<p/KTe 

=0 r-2 = ® - "i) (■'•9) 

o X 

where 

Sq = permittivity of vacuum, 
n^ = unperturbed electron density 
K = Boltzmann constant 
Te = electron temperature, and 

other variables are as in Eqs,(l,7) and (1.8). 


The following dimensionless variables make all 


the coefficients unity; 
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where 


= xAd 

= Op 

= e9/K Te 




e f. K T 

Debye length = — 9 — ) ' 


plasma ion frequency 


acoustic speed (T^ = 0) 


M ^ 




ion mass 


The Eqs.(1,7) to ( 1 ,9) now become; 


4- V ' 


( 1 . 10 ) 




(n’v' ) 


( 1 . 11 ) 


a 


a ~ n’ 


(1 . 12 ) 


Now, n' , X and v’ can be expanded in terms of a 
dimensionless amplitude parameter : 


1 + 6n^ +6 n^ -s- 


6X^ + 6^ X^ + ... 


(1 .13) 


6v^ + 6 Vrt + ... 
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We also transform the variables x’ and t' to li and % 
through 5 


. 1/2 


(x* - t* ) 


= t' 


(1 .14) 


So that 


8 

8t' 


53/2 




8 

8T 


8 


6 V2 


8 

51 


(1 .15) 


If we now assume that n^ ^ I -j » tend to zero as 

I § 1 -* oc t it can be shown with the help of Eqs,(1 ,10) to 
(1 .15) tha t [ 5 ] s 


n 


1 



= U 


'Atiere U satisfies the KdV equations 


au 

It 


+ U 


3 U 


1 

2 


3^U 

8g3 


0 


(1 .16) 


(1 .17) 


Thus ion waves of amplitude one order higher than linear 
are described by the KdV equation. 


.4 Literature Review 

The soli ton was discovered and named in 1965 by 
Vjabuskey and Krauskal [ sl , v\ho were experimenting with the 
imerical solution of the KdV equation. The KdV equation 


I 
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had been introduced earlier at the end of the last century 

to describe wave motion in shallow canals by D,J. Korteweg 

and G, de Vries (see Section 1*3,1), The first scientific 

description of soliton as a natural phenomenon goes however 

of 

back to first half/ the nineteenth century and was reported 
by J. Scott-Russel in 1844 [4,9,10,1l]. 

But the real breakthrough occurred in 1967 when the 
idea of spectral transform technique was introduced by 
Gardener, Green, Kruskal and Miura as a means to solve the 
Cauchy problem for the KdV equation [12] , Soon afterwards 
Lax put the method into a framework, that provided a clear 
indication of its genera lity and greatly influenced future 
developments. A few years later, in 1971, Zakharov and 
Shabat, by a nontrivial extension of previous approaches 
solved the Cauchy problem for another important nonlinear 
evolution equation, so-called nonlinear Schroedinger 
equation [9,13] . The way v£is thereby opened for the 

search and discovery of several other nonlinear evolution 
equations, or rather classes of such equations, solvable by 
these techniques. This process continues unabated to the 
present day. The applications of this subject have been 
percolating through the whole of physics (from nonlinear 
optics to hydrodynamics, from plasma to elementary particle 
physics, from lattice dynamics to electrical networks, to 
superconductivity, to cosmology and to epidermiology etc). 



This is of course related to the central role played 
by non-linear evolution in mathematical physics, and • 
more generally in applied mathematics. In fact the spectral 
transform approach (see Section 2,2) constitute in some sense 
an extension to a nonlinear context of the Fourier transform 
technique, vtiose all-pervading role for solving and investi- 
gating linear phenomena is well known. 

We will now briefly outline the numerical simulation 
studies for nonlinear evolution equations reported in 
literature. As has been mentioned above Zabusky and 
Kruskal [8] were the first to study KdV equation through 
a finite difference method. They used a leapfrog scheme 
(see Section 3.2). An improved scheme was later obtained by 
Grieg and Morris [14] , This was a hopscotch scheme, 

Fornberg and Witham [25] devised an accurate scheme using 
spectral methods for x-variable and leapfrog in t. 

Among the finite element methods, Wahlbin [15] suggested 
a dissipative Galorkin scheme. Other workers v\ho used 
Galerkin or Petro-Galerkin schemes are Andersen and Michell [26], 
Alexander and Morris [16] , Winthev, Mitchell and Schoombie 
[17] and Sanz-Serna and Chriesti [18]. 

1 .5 Outline of the Present Work 

In Chapter 2 we describe some analytical methods to 
obtain soliton solutions of the KdV equation. In Section 2,1 
a method to construct single-soliton solutions is described, 
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whereas in Section 2,2 a more general technique (the 
Spectral Transform Method) is briefly outlined. 

In Chapter 3 the two main numerical methods, viz, 
the finite difference and the finite element techniques 
are briefly outlined for partial differential equations 
in general and for the KdV equation in particular. This 
is followed by the actual numerical computations and their 
results in Chapter 4, The results obtained during this work 
are compared in detail with the previously available 
results in Chapter 5, 
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CHAPT ER 2 

ANALYTICAL METHODS FOR THE KdV EQUATION 

Sinqle-Solit on Solutions for the KdV Equation 

In order to obtain single-soliton solutions of the 
KdV equation we seek travelling solution of permanent 
shape and size by trying solutions of the form 

u(x,t) = U( 5), (2.1) 

v\tiere 

^ = X - ct, ( 2, 2) 

As a result we have, 

Ut = -c U§ 

= Ug (2.3) 

'^xxx ^ 

Substituting Eqs,(2,1) to (2,3) in Eq.( 1 .3) we get the 
ordinary differential equation 

-c + U Ug + e 0, (2.4) 

which can be integrated once with respect to ^ to yield 

-c U + + e Ugg = A , (2,5) 

where A is a constant of integration. Multiplying Eq.(2,5) 
by U g and integrating we get. 
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- 5 c + I + 2 "" AU + B (2.6) 

For soliton solutions we use the boundary conditions, 

U, Ug, Ug^-* 0 as ISI -* 00 . Therefore, we have 

A = B = 0 

With these values of A and B and with a slight rearrangement, 
Eq.( 2.6) can be written as 

2 

Ug = ^ .(3c - U) (2.7) 


taking square roots of both sides in Eq,(2.7), the dependent 
variables U and the independent variable t may be separated 
as follows: 


Nf . 3. e dU, 

U -“u 


( 2 . 8 ) 


The last equation may be integrated using standard 
methods to yield 


U(g) = 3c Sech^ (^ d), (2.9) 


where c is a constant > 0 and d is an arbitrary constant. 
By substituting for ^ from Eq,('^.2) we obtain 

u(x,t) = 3c Sech^ (2,10) 

This is a single-soliton solution of the KdV equation,iMl 7 rch 
we have been seeking. It can be seen that the amplitude 
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width and speed are all determined by a single parameter c, 

and are thus interdependent. In particular the pulse height, 

-1 /2 

width and speed are proportional to c, c , and c respecti- 
vely. 

Lastly we may remark that if we put e = 1 in Eg, (2,10) 
we get the single-soliton solution for the KdV Eg. (1,4): 

u(x,t) = 3c Sech^ (^ nTc (X - ct) + d) (2.11) 

Furthermore, in Egs.(2,10) and (2,11j if we take d = 0, the 
peak of the soli ton coincides with x = 0 at t = 0, A plot of 
Eg, (2.10) for d = 0 and c = 1 at t = 0 is shown in Fig. (2.1) 
for two different values of e . 

2,2 Spectral Transform Method for the KdV Equation 

In this section we briefly discuss Schroedinger spectral 
transform which is appropriate for solving a class of nonlinear 
evolution eguations of which KdV eguation is a special case. 

For this purpose we consider the stationary Schroedinger 
eguation: 

u(x)’i'(x,k) = k^^(x,k) -oo<x<+ oo 

( 2 . 11 ) 

Here u(x) is a gi^^’en real function that vanishes sufficiently 
fast asymptotical||y, say 




( o‘ X) n 


FIG. 2.1 A PLOT OF A KdV 50LIT0N Eq-2.10 WITH d=0 ATT 
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Eq,(2,11) is a singular Sturm-Liouville problem. The spectrum 
of the Sturm-Liouville problem characterized by Eq,(2,11) 
consists of two components; 

1 , A continuum, including all positive values of 

2 

the eigen value, k , 

2. A number of discrete negative eigen values 

k = ~ Pn» Pn = "I »2, . . . .,N. 

The continuum part of the spectrum is completely characterized 
by the reflection coefficients R( k) [9,193 . On the other hand 
to completely characterize the discrete part one needs not 
only the values of but also the corresponding normalization 
coefficients given by; 

-!-00 O —1 

= [ / dx fn^(2:)3 , n = 1,2,...,N 

-00 

(2.13) 

v\hGre 

f^(x) is the eigen function corresponding to p^. 

The spectral transform S of the function \i(x) is by 
definition the collection of data 


S( u) = { R( k ) . - 00 < k < + 00 ; P » n = 1 , 2, . , . . , N } 

^ n 

( 2 . 14 ) 

one-to-one 

It is found that there is £ correspondence between functions 
u(x) and the spectral transform S [u]. If u also depends on 
another variable say t, i.e., u =u(x,t), then we get a 


time dependent spectral transform; 

u( x,t) s( t) 


( 2 . 15 ) 



The all-important, and highly nontrivial discovery of 
the last two decades has been that there exists a class of 
physically interesting time evolutions of u(x,t) governed by 
nonlinear partial differential equations for which the time 
evolution of S( t) is governed by linear ordinary differential 
equations. Thus the time evolution is much simpler in the 
spectral space than in the configuration space. Once the 
problem has been solved in spectral space the inverse spectral 
transform can be obtained by Vw'ell defined procedures, yielding 
the final solution. As is clear from this discussion this 
procedure is analogous to Fourier technique for solving certain 
linear problems. However, the details mil not be pursued 
any further here [9,20,21]. 

In the next Chapter vje will discuss the methods for 
obtaining the approximate numerical solutions to the KdV 
equation. 
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CHAPTER 3 

NUMERICAL METHODS APPLICABLE TO THE KdV EQUATION 

3.1 Introduc ti on 

As mentioned in earlier chapters, the KdV equation is 
a nonlinear partial differential equation having localized 
travelling solutions, called solitons. All nonlinear wave 
equations with soliton solutions are solvable by the powerful 
tool of Spectral Transform as discussed briefly in Section 2.2. 
There are many other problems, hovjever, of equal interest for 
which the Spectral Transform method cannot be applied . Indeed, 
it is in this area that many of the unsolved problems related 
to the soliton - like behaviour occur. 

Therefore, a useful approach vdll be to study the 
evolution of solutions numerically. Numerical methods can 
be tested on equations for which analytic results and exact 
solutions are known, and then applied to equations for which 
analytical results are not known. 

3.2 Basic Numerical Metho ds 

As a tool for the study of partial differential 
equations, numerical analysis offers a confusing variety of 
technir ues. For a specific exfuation, the question of the 
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required for a specific range of the independent variables, 
the limitations of time and storage of space and machine 
word length. Equally important is the amount of time and 
effort available for the development of the appropriate ■ 
s oftvjare , 

Most of the nonlinear evolution equations can be put 
in the form 

= L(u), or (3,1) 

= L(u) (3,2) 

vhere L(u) is some general non-linear differential operator in 
space variable x. 

In order to treat the equation numerically we must 
replace boundary conditions ’at infinity by conditions at some 
finite boundary. To minimize the effect of an infinite boundary, 
we can take the boundary at large distance from the regions 
where u is nonconstant or, where it effectively vanishes . At 
these boundaries we take u and its space derivatives as zero. 

Now the first step in solving a differential equation is 
to discretize it, i.e,, the differential equation must be 
replaced by an approximately finite system of algebraic equations. 
There are three basic discretization methods; 

1 , Finite Difference methods 

2. Finite element methods 


3, The Method of lines. 

Here we will discuss the first two only because of their 
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relevance to the present, problem. 


3 .3 The Fi nite Differen c e Methods 

The main idea behind the finite difference method for 
obtaining the solution of a given differential equation is to 
approximate the derivatives appearing in the equation by a 
set of values of function at selected number of points 
called nodes. The most popular way to generate these approxi- 
mations is through use of Taylor series, but other approaches 
can be used. One possibility is the direct application of 
the interpolation formulae. 

If Xj_ + h = » then from Taylor series, we have 

Df(xjj_) = + 0(h) (3.3) 

v\here 

^ -- dx * 

Another finite difference approximation is the central 

difference as opposed to the forward difference approximation 

given by Eg. (3. 3), Given below is a central difference 

2 

approximation to the second order derivative, D f: 


D^f(Xj_) 


f(xi.,^) - 2f(xi) + 

^ 

h 


(3.4) 


Using either the Taylor series or the interpolation 
formulae, the finite difference approximations to various 
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derivatives of functions of one or more variables can be 
constructed. The use of these approximations in the 
differential equation together with the initial boundary 
conditions yields an algebraic system of equations which can 
be solved, using aporopriate techniques. Such schemes have 
been used for the nonlinear KdV equation by Zabusky and 
Kruskal [ s] , and Grieg and Morris [14], 5A?hereas, the results 
obtained in this work (see Chapter 5) are compared with their 
work, we will, however, not discuss the details of their 
schemes . 

3 *4 Finite Element Metho ds 

Let the partial differential equation be of the form 

Lu = f (3.5) 

where L is the differential operator, then the finite element 
method involves the follovang steps; 

1 , Select some known basis functions to approximate, 

say u(x) in one dimensions 

N 

u(x) ^ s Qtj 9 j(x) = u(x,a) (3.6) 

3=1 

Here 9j(x) are basis functions, and a denotes the set of 
a i*s , 

2. Substitute U(x,a) into Eq,( 3.5)5 then the 
unknown coefficients will be approximated by solving 
the system LU = f. 
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3. Since stop 2 is generally impossible to do 
exactly (one has only n unkno\wis available to satisfy equations 
viiich hold for all values of x in some domain), select N 
suitable conditions to obtain the N unknowns oc^. 

The method outlined above is more general than the 
finite element method. The finite element method requires, 
in addition that the basis functions bo finite elements? i.e., 
functions which are zero except on a small part of the domain 
under consideration. The piecewise polynomials are typical 
examples of finite element basis functions. The most commonly 
used basis elements are piecewise linear. 

To obtain the unknowns a^’s, one proceeds as follows: 
The approximate solution U = S Eq„(3,6) is 

substituted in Eq,(3,5), Since U is not an exact solution, 
there results a residual 

R(x,a) = LU - f (3.7) 

which is generally not equal to zero. Our aim will be to 
select the coefficients oc^ in such a way that the residual 
is forced to be zero in some average sense. This can be 
accomplished by selecting a set of weighting functions w^’s 
and then setting the integrated weighted residuals to zero. 

We write it s'/mbolically as 

(R, w^) = /^R(x,a) Wj dQ = 0, (j = 1,2, ....N) (3.8) 


Here Q is the domain of interest. 
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The solution of Eq*(3,8) vdll provide us the values 
for the set of M unknown coefficients. The weighting function 
Wj can be chosen in many ways^- Following approaches are 
generally adopted, 

3.4.1 T he Collocation Method 

An obvious choice to determine a^^'s is to make 
LU( ) = 0 at N selected points. Here the weighting functions, 
Wj’s, are taken to bo shifted Dirac Delta functions, i,e,, 

Wj(x) = 6(x~Xj), j = 1,2,....,N (3.9) 

This function 6 ( x-x j ) has the important property that 

f( x) 6 ( x-Xj )dx = f(xj) (3.10) 

cond the coefficients are determined from the conditions 

/^R(x,a ) 6(x-x,-) dQ =0, j=1,2, (3.11) 

The main disadvantage of the method is that the accuracy 
obtainable depends strongly on the location of [22,23] , 

3.4.2 The Least Square Method 

In this method the residual itself is taken as the 
weighting function, i,e,, 

Wj = R(x,a) . (3,12) 

Then the coefficients are determined by minimizing the integral 
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1 = Sq R^(x,a) dQ = /q (LU - f)^ dQ , (3,13) 

that is, 

=0, j = . (3.14) 

j 

Here the integral of square of the error is minimized, 
hence the name, least square method, 

3 • 4 .3 The Galerkin Method 

Perhaps the best known method based on the weighted 
residual principle is the Galerkin method. Here the weighting 
function is chosen to.be one of the basis function itself, 
i.e, w.(x) = 9 -(x), and the constraint imposed for the 

sj sj 

determination of the coefficients ® j ' s is s 

(R(x,a ) » 9j) = 0 

or 

(LU - f, 93 -) = 0 (3.15) 

Since 9 j's (basis functions) form a complete set, the Galerkin 
method can be interpreted as forcing all residuals to be 
orthogonal to the specified basis functions. Therefore good 
numerical results can be expected and are in fact obtained. 

3 .5 A Comparison of Different Methods 

Now a choice of method will depend mainly on the 
nature of the problem to be solved, i.e. the order of the 
differential equation, boundary conditions and continuity 
irequirements , If one has a rectangular domain, then the finite 
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element methods becomes simpler and more efficient to set up. 

If one has a more complicated domain, the finite difference 
method becomes more complicated and less accurate, v\tiereas 
the finite element methods are relatively unaffected. The 
principal advantage of the finite element methods is that they 
are generally more accurate. The second advantage, as indicated 
above, is that they are more flexible in handling complex 
geometrical shapes by the use of arbitrary shaped simple 
elements , 

3*6 C onvergence, Consistency and St a b ility 

V\ihatGver approach is chosen, one would ideally like to 
show that the method is convergent [24], i,e, the global 
error for some fixed and finite value of time tends to zero 
as the step length (nodal separations) tend to zero, or as 
the number of basis functions tend to infinity. Unfortunately, 
the proof of convergence is usually extremely difficult, even 
when possible, and a satisfactory theory is only available in 
general for linear equations , Most discussions of nonlinear 
methods are based on an analysis of the equation linearized 
about some constant solution. 

To ensure convergence, methods in general must be 
consistent and stable. A consistent method is one in vhich 
the truncation error tends to zero as the step lengths in 
the problem tend to zero. Sometimes this limit must be taken 
in a specific way; for example in the Hopscotch method the limit 
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must be such that k/h -» 0 as k, h -♦ 0, v>/here h and k are 
step lengths in space and time respectively. Stability means 
that some norm of the approximate solution must remain bounded 

as n CO , where nk = t for a fixed t. Even stability is 
usually difficult to prove for nonlinear equations, and the 
usual approach is to study the stability of a linearized 
version of the equation [16,24]. 

In the next chapter, we discuss the development of a 
finite element computer code for the KdV equation using cubic 
splines [27,28] . As will be seen there, the use of finite 
elements in space for the KdV equation yields a nonlinear 
system of ordinary differential equations which are solved by 
a finite difference method. The results obtained will be 
presented and discussed in latter chapters. 
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CHAPTER 4 

NUMERICAL CO?APUTATIONS AND RESULTS 

4.1 Preliminary Remarks 

In our work we have used Finite Element Method (FEM), 
and in particular the Galerkin method using cubic splines 
for the solution of the KdV equation. The KdV equation is 
generally written in the forms 

[u] = ^ ^xxx ° 

where e is a positive constant (see Section 1,2), As 
mentioned in earlier chapters it possesses the single-soli ton 
s olutions 

u(x,t) = 3c Sech^ ( kx - wt+ d) (4,2) 

where 

k = ^ (g')"'^^, and 

w = kc (4.3) 

Here c and d are parameters, c is the velocity of the soliton, 
and d is the phase or central shift, 

4.2 Calculation of Galerkin Constants 

We take the interval of study as [0,2] on the x-axis, 
because at the boundaries of this interval the solution almost 
vanishes for the values of c, and d for which we will be 
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interested and over the time range of interest. 

Now the interval [ 0j2J of the x-axis is divided into 

o 

N subintervals or elements each of v^/hich has a length h = , 

We approximate the solution as 

N-2 

U = S a (t) <P.(x) (4.4) 

j-2 ^ ^ 

viiere the basis functions cpj(x) are approximately shifted 
basic cubic splines. These will have the local support 
on the interval [(j-2)h, (j+2)h]. 

Now substituting (4.4) into (4.1), we get the residual 

R(x,a) = + UU^ + (4.5) 

As mentioned in Chapter 3, in the Galerkin method, we take 
the Vifeighting functions also as 9j(x) and impose the 
condition: 

2 

/ R(x,a) 9-;(x) dx = 0, j=2,3, . . . ,N-2 

Q 

(4.5) 

where the detailed expressions for R(x,oc) is seen to be 

^(x;,a) = ^ ( ? aj(t),9j(x)) + ^ “j^t) ‘Pj(x). 

|j( Saj(t).9j(x)) + 

(4.7) 

where the summation over j is fron 2 to (N-2) as in 
Eqs .(4.4) and (4.5) » 
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Substituting Eq,(4,7) into Eq,(4,6) vve get the 
following system of coupled ordinary differential equations 
for the time-dependent expansion coefficients, a^lt) 

ba .. . 

“j“k - h-tc..)a. =0 (4.8) 

Here the summation over the repeated incjices (i.e. if they 
are repeated in the same term) is implied from 2 to (N-2), 

The Galerkin constants A. B. and C . • are given by 

IJ IJK IJ 

+ 00 

^ij ^ 9i(x).9j(jc) bx 

+00 

%jk " (9i»9j9^i,^b = ^ 9j_(x).9j(x) 9^i^^\x) dx 

+ 00 

Ci. = h^( ^)9 ^.^-0 = h^ 1 9i(x) 9^-^\x)dx 

J ^ -00 ^ 

= 9^.^\x)9/^\x) dx 

-CX) ^ 

(4.9) 

Since our region of interest is [0j2], the above integrations 
are done in the interval [0,2], Furthermore, the basic 
cubic spline functions 9 j(x) are spread over a compact local 
support of 4h only. Hence the integration exists over an 
interval of 4h giving the band nature of the matrix. 


31 


As mentioned in the beginning, the basis functions 

9 j’s have the support over the interval [(j- 2 )h, (j+ 2 )h] , 

Therefore, 9 ^ vail be spread over the interval [-h, 3h] , 

we 

Since we have taken our interval as [0,2] ,/take 92 as the 
first basis function. Similarly our last basis function vail 
be 9 j ^_2 and, not t because the latter will extend to one 

step beyond the interval, 

Nov/, we write a basic cubic spline in the following 
standard form [16,2^ 

M( y) = 0 y < “2 

= l[(y+2)^] -2 ^ y ^ -1 

= ■5[(y+2)^-4( y+i )^] “1 1 y 1 0 ( 4 . 10 ) 

= I [(-y-i-2)3-4(-y-i-l)^] 0 < y < 1 

= -5 [(“-y+2)^] 1 < y i 2 

= 0 2 < y 

Then the basis functions can be written as 

9j(x) = M(| - j) (4.11) 

Once the basis functions are known we can integrate (4,9) to 
get the Galerkin constants. 

In calculating these constants, the first thing. 

We will need is a subroutine to multiply two polynomials 
because of the product term present in the integrand. The 
product polynomial can then be integrated easily by suitably 
changing the power and the coefficient of each term automatically 
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by the computer code. The integration can further be simpli- 

of 

fied by changing the limits/ integrations to [0,1J in all 
cases by suitable rescaling of the variable of integration. 
Since the polynomials are defined piecewise on intervals 
of h, the integrals can always be reduced to the forms 

1 P 

/ [ a^ + a^y + a y’^] dy = 2 k +1 ) 

o ^ k =0 

(4.12) 


4 *3 Initi al Conditions 

After the calculation of the Galerkin constants, we 
will need a subroutine to solve the system of coupled 
nonlinear ordinary differential equations, (4.8). Hence 
the initial value of the dependent variables, a j( O) , need to 
be specified. Here we will use least square technique to: 
calculate the initial values. We minimize the expression 

2 

^ f; I aj(0) <pj(x) - ^^(x)!!^ 

2 

and 

= ( ? " ^o’ k “ k^^^^k^ 

(4.13) 

where f^ = u(x,0), and u(x, t) is given by Eq,(4,2). The 
particular numerical values of , c and d needed to evaluate 

u(X{,0) are given below, 

3 E 

Setting — = 0, we obtain the following equations from 

aeCjAu; 

which to compute a j(0) 2 
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5 

3 





9i)/h 


(4.14) 


where can be seen to be the same as given by the first 

of Eq.( 4,9) . 

Since the integral on the right hand side of Eq,(4,14) 
cannot be easily evaluated analytically, we have used a 
numerical integration to perform the required integration. 


Here the matrix formed by A. . is a band symmetric 

J 

matrix of order (N-3) . We can save much space in storage by 
utilizing its band nature, i.e, storing it in band symmetric 
storage mode. The number of co-diagonals of the symmetric 


matrix formed by A^j is 3 here 


4.4 Solution of the ODE System 

The Eqs,(4.S) may be written in the form 


A 


ij 


da^. 

dt 


h-3 e.jaj - j ^ (4.15) 


which is equivalent to the following matrix form; 


A -gY = F(a,T) (4.16) 

where A is the band symmetric matrix and F(°^,T) is a column 
vector. In the present case, however, F doos not depend 
explicitly on time. Since available differential equation 
solvers [IMSL, NAG Libraries] , generally requires the system 
of equations in the form = F(a,T) , we need to convert 
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Eq. (4,15) in the above form by calculating A and then 
writing in the desired form; 

Da 

UT ~ ^ F(a,T) (4,17) 

The system of Hqs,(4,17) is then solved, using standard 
subroutines LIN2PB and DREBS from IMSL library. 


^ ^ Evaluation of the Cub ic Spli ne for the Final Solution 
Once the a (t), j = 2, ,N~2 have been obtained 

j 

at the desired times T by solving the system of equations 
(4,17), the values of U(X,T) can be evaluated using Eqs,(4,4) 
and (4.1 1 ) . 

Figure (4.1) gives an overall outline of the computer 
code developed for carrying out the above-mentioned tasks. 

The details of the subroutines used are given in Appendix B, 
and the results of some sample calculations are presented in 
Appendix C, In the following section, we give a summary of 
these calculations, particularly with a view to bring out 
the accuracy of the computed solutions, 

4,6 Results and Their Accuracy 

The actual numerical computations were done for the 
following values of the parameters in Eq. (4,2) ; 
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e = 4.84x10 

c = 0.3 

d = ^ and 

w = kc , 


(4.18) 


These values were chosen to facilitate comparison with 
previous works. For the same reasons the interval of 
interest was taken as [o,2]. It can be seen from the actual 
initial data as well as its cubic spline representation 
given in Table (4.1), that at the ends of the interval the 
initial values are negligibly small and therefore can be 
taken as zero. Furthermore we will restrict the evolution of 
the solution to a maximum of T = 2.0 second. Upto this time, 
it can be easily seen that the single-soli ton solution 
continues to remain negligibly small at the end points. 

We have calculated the solution for two different 
grids, N = 40 and N = 50, In the former case the final 
solution was computed for T = 0,5, 1,0, 1,5 and 2.0 seconds. 
For N = 50 tho solution was obtained for T = 0.4, 0,8 and 
1,2 second. The results are summarized in Tables (4.1) to 
(4,6), and Figs. (4.2) and (4.3). 

Table (4.1) shows at alternative grid points the 
values of the cubic spline approximation ' of the initial data 
for N = 40, It is seen that an maximum error of approximately 
0,4 percent is incurred v\hen compared to the peak value. 
Therefore, one perhaps cannot expect an accuracy better than 
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this in the final solutions, (This is confirmed later). 

The corresponding error for N = 50 is seen to be about '' 

0,2 percent. 

In Fig. (4,2) the initial data at T = 0,0, the calcu- 
lated and expected solutions for T = 2,0 are plotted. The 
figure clearly brings out the overall soundness of the 
computer code developed. It is infact difficult to distin- 
guish between the calculated and expected solutions graphi- 
cally, The same is true for Fig, (4,3) where the expected 
and calculated solutions are shown for T = 1 ,2 using M = 50, 

Tables (4.2) and (4,3) show the calculated and 
expected positions of the single-soliton peaks at different 
times for N = 40 and 50 respectively. It is seen that in 
all cases the peaks are observed precisely where expected,. 

In Table (4,4) a detailed comparison of the calculated 
and expected solutions, and the corresponding absolute and 
relative errors at alternative grid points, are given. The 
maximum error expressed as a percentage of the peak value 
is seen to bo approximately 2 percent. Here it may be recalled 
that the corresponding error in the initial data represen- 
tation was about 0.2 percent (see Table 4,1 and the discussion 
above). The and L 2 errors are also given at the bottom 
of Table (4,4), 

Tables (4,5) and (4,6) give the errors in the peak 
values only for various runs. It is seen that the error 



is generally less for N = 50 than it is f or N = 40, Similarly 
it generally increases if time of evolution is increased. It 
should be pointed out, however, that it is and L 2 errors 

which give a more systematic idea of the accuracy of calcu- 
lations. These will be preserited in Chapter 5, v\here a 
detailed comparison will also be made with previous works. 
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AN OUTLINE OF THE PRESENT COMPUTER CODE 
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TABLE 4,1 

Comparison of Actual Initial Data and Its Cubic Spline 

, . * 

Representation (N = 40) 


s. 

No. 

X 

Actual Initial Data 

Uq(X) Cubic Spline 

Representation U^CX) 

1 

0.20 

0.80608E-08 

-0.2381 3E-06 

2 

0.40 

0.1 171 8E~05 

-0.19135E-05 

3 

0.60 

0.17033E~03 

0.13309E-03 

4 

0.80 

0.24427E-01 

0.24024E-01 

5 

1 .00 

0.90000E+00 

0.90392E+00 

6 

1 .20 

0.24427E-01 

0.24024E-01 

7 

1 .40 

0,17033E-03 

0.13309E-03 

8 

1 .60 

0.1 1718E-05 

-0.19135E-05 

9 

1 .80 

0,80608E-08 

-0.23808E-06 

* 

Notes 

s 1 . 

error (N=40) = 3 

.92E-3, error (N==40) =1 .77E-: 


2. 

L and Lo errors in 

found to be 2,04E-3 ■ 

N=50 representation were 
and 0,81 E~3 respectively. 



TABLE 4.2 

Ejrrors in th© Observod Locations of th© Soliton Peaks 
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Errcrs in Peak Values for Different Evolution Times (T) 
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Errors in Peak Values for Different Evolution Times (T) 
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CHAPTER 5 

COMPARISON WITH PREVIOUS WORKS AND CONCLUSIONS 

5*1 Comparison of L^and L 2 Errors 

A variety of finite element schemes and some 
finite difference schemes have been used for the numerical 
solution of KdV equation by previous workers . Here we 
compare our results with some of the previously available 
results , 

Alexander and Morris [16] have used a Galerkin 
scheme using cubic splines for basis functions and an 
appropriate set of quintic boundary functions. They have 
also used a dissipative Galerkin scheme and presented their 
results for different values of dissipation parameter, q. 
However, they have evolved the single-soli ton solution 
upto 0,4732 second, v\hereas, we have evolved the solution 

uPto 2 seconds . Furthermore, they have given only L errors 

00 

in their calculations. In Table (5,1) we compare our results 

with the above-mentioned calculations. We find that the 

L error in our run for T = 0.5 is less than the error for 
00 

q = 0 and T - 0,3958 in [ 1 6] , The errors in other runs 
also compare somewhat favorably for the present work. 

In Table (5,2.) we compare our calculations with 
N = 50 vdth those of Alexander-Morris [16] with N = 60. 

Here the results of the present work appear to be appreciably 
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better as can be seen from the table. In particular when 

we evolve the solution to T = 1 .2, an L - errors of 

00 

-3 . . 

17,3x10 is incurred v\hich is almost comparable to the 
Loo“en‘or in [l6], v\tien the solution is evolved only to 
0,4617 second. 

In Table (h,3), the errors are compared with the 
finite difference schemes of Zabusky-K3?uskal [8] and Grieg- 
Morris [14], Here both the and errors in the present 
work are found to be a factor of 2 to 10 times less in the 
present work. 

Lastly in Table (5,4) we compare our results with 
some recently used Petrov-Galerkin schemes [17], Here we 
find that the - errors in the present calculations are 
about a factor of 4 less than those in the Petrov-Galerkin 
schemes. However, one of the later schemes appears to have 
somev\tiat less L^^-errors , We also present a rough comparison 
of computation times in Table (5,5), 

5,2 Co nclusions and Suggestions 

In conclusion we may say that the present attempt to 
develop a computer code for the numerical evolution of 
nonlinear wave equations has been highly successful and the 
results obtained already compare quite favorably with the 
previous works. However the present code needs to be devel- 
oped and generalized in many directions, or a fresh more 
general code can be written based on the experience obtained 
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here, we list some of the points below; 

1. The exponential, trigonometric or hyperbolic 
splines may be used. An exponential/hyperbolic polynomial 
may be particularly useful in the present context as a 
boundary function. 

2. A variable/adjustable grid can be used, 

3. The code should be written such that it can 
handle a variety of basis and weighting functions, as well 
as assemble the Galerkin constants for many nonlinear 
evolution equations of interest, 

4. The finite element method can be used for time 


variable also 
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APPENDIX A 


TuC TR AN3FQE ivu\TI0NS FOR £Q. (1.1) TO EQ. ( 1 .4) 


iy trying the transformations 


u -* a u, 

t {3 t, and 

X -* Y 

Viherv:' a,p and y are constants, one finds that the follov\/ing 
V’lu.r; of a, p and y transform Eq.( 1 .1 ) toEq.(1,4)s 
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